function [Grid,TrMat]=AR1Approx_Rouwenhorst(N,Ybar,Rho,Sig)
%% Preliminary
% CHECK IF abs(rho)<1 
if abs(Rho)>=1
    error('The persistence parameter Rho is not smaller than 1!');
end

% CHECK IF n_R IS AN INTEGER GREATER THAN ONE.
if N<1.50001
    error('The number of grid points N must be greater than one.');  
end

% CHECK IF n_R IS AN INTEGER.
if mod(N,1)~=0 
    warning('The number of grid points is not an integer! Rounded value will be used!')
    N=round(N);
end

%% Construct the Grid and Transition Matrix
p = (1+Rho)/2;
q = (1+Rho)/2;
psi = sqrt(N-1)*Sig/sqrt(1-Rho^2);
mu = Ybar/(1-Rho);

Grid = mu+psi*linspace(-1,1,N)';
TrMat = [p 1-p;1-q,q];
ii = 2;
while ii<N
    Mat_1 = [TrMat,zeros(ii,1);zeros(1,ii+1)];
    Mat_2 = [zeros(ii,1),TrMat;zeros(1,ii+1)];
    Mat_3 = [zeros(1,ii+1);TrMat,zeros(ii,1)];
    Mat_4 = [zeros(1,ii+1);zeros(ii,1),TrMat];
    TrMat = p*Mat_1+(1-p)*Mat_2+(1-q)*Mat_3+q*Mat_4;
    TrMat(2:end-1,:) = TrMat(2:end-1,:)/2;
    ii = ii+1;
end
TrMat = TrMat';
TrMat = bsxfun(@rdivide,TrMat,sum(TrMat));
